Screening of charged spheroidal colloidal particles 
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C . 

r-{ , We study the efTective screened electrostatic potential created by a spheroidal colloidal particle 

immersed in an electrolyte, within the mean field approximation, using Poisson-Botzmann equation 
in its linear and nonlinear forms, and also beyond the mean field by means of Monte Carlo computer 
simulation. The anisotropic shape of the particle has a strong effect on the screened potential, even 
at large distances (compared to the Debye length) from it. To quantify this anisotropy effect, we 
f~| ' focus our study on the dependence of the potential on the position of the observation point with 

O ' respect with the orientation of the spheroidal particle. For several different boundary conditions 

, (constant potential, or constant surface charge) we find that, at large distance, the potential is 

S ■ higher in the direction of the large axis of the spheroidal particle. 
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INTRODUCTION 



Charged colloidal systems are widespread in nature and their study is important as they include systems such as 
'"^ ■ clays and protein solutions in water. Derjaguin, Landau, Verwey and Overbeek studied the stability of lyophobic 
^ ' colloids by solving the Poisson-Boltzmann equation [1, 2] and stated the theory known as DLVO. They found the 
O . screened electrostatic potential for a spherical colloidal particle of radius R and valence Z, immersed in an electrolyte 
with positive ions of charge z+e and density n+, and negative ions of charge — z_e and density n_. In the linearized 
(Debye-Hiickel) regime, the screened potential is given by 



^ ' pK,R p-K,r 

00 

where ^ denotes the dimensionless potential /3ex potential, Ib = /Se'^/e is the Bjerrum length, e is the elementary 
charge, /3 — {kBT)~^, with T the temperature, ks the Boltzmann constant, and = AirlBiz^n^ + z^n^) is the 
inverse Debye length. 

However, many colloidal systems of importance have, in their disperse phase, macromolecules whose shape cannot 
be considered spherical, for example, mineral liquid crystals [3], and thus their behavior is not properly explained 
by the solution of Poisson-Boltzmann equation for the case of a spherical coarse grained particle. Therefore it is 
worthwhile to study models for non spherical particles, for example ellipsoids of revolution (also known as spheroids). 
Furthermore, there are now synthetic clays with uniform particle shapes like laponite [4], which can be modeled as 
\^ ' oblate spheroids and it is also possible to synthesize spheroidal nano-particles [5], which provides experimental data 
5^ , for these models. 

The solution of the linear Poisson-Boltzmann equation has been carried out for particles with shapes other than 
spherical (or the infinite plane case), in particular the cylinder [6, 7], disc [8-10] and also the more general case of 
ellipsoids of revolution [11, 12] have been studied. For particles with a disc shape, the screened potential has a long 
range behavior which decays with the distance as the Yukawa potential, but multiplied by an anisotropy function [8], 
which characterizes the angular dependence of the electrostatic potential at large distance. This is actually a general 
feature for anisotropic charged objects immersed in an electrolyte: the screened electrostatic potential is anisotropic 
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even at large distances from the charged object [13]. This can be contrasted to the situation of a charged object in 
the vacuum, where the electrostatic potential created by this object becomes isotropic at large distances: the first 
term in the multipolar expansion is isotropic. In this work, we study further this anisotropy function for the general 
case of charged ellipsoids of revolution in the Debye-Hiickel as well as in the nonlinear Poisson-Boltzmann regimes. 
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FIG. 1: Prolate (left) and oblate (right) spheroids. 

For the present geometry of the system the spheroidal coordinates will be used. The spheroidal coordinates 77, (j)) 
are related to the cartesian coordinates {x,y,z) by [14]: 



X 


= aVie 


-1)(1 




y 


= aVie 


-1)(1 




z 








X 






- T? ) 


y 


= aVie 






z 









COS( 



(2) 

for the prolate case, and 



(3) 

for the oblate case. Here 2a is the distance between the foci of the ellipse. The coordinate ^ is called the radial 

coordinate since it plays a similar role as r in spherical coordinates (r, 0,(j)). The coordinate rj will be called angular 
coordinate, as it plays a similar role to cos 9 in the spherical case. The surface of constant ^ is a spheroid, with major 
and minor semi-axes, R and R' given by 

R = a^, 



R' = av/^2T:i (4) 

for the prolate case, and 

R' = a^, 

R = + 1 (5) 

for the oblate case, see figure 1. Note that the change of variables ^ — >■ ^-nd a — > —ia in equations (2), transforms 
from the prolate to the oblate coordinates, so in the following we will treat in detail the prolate case, and only 
report the results for the oblate case, without repeating the calculations. In section II, we study the solution to 
Poisson-Boltzmann equation in the linearized regime, for spheroidal colloidal particles in two particular solvable 
cases: constant surface charge for an ion-penetrable colloid and constant surface potential. Section III deals with 
the numerical solution of the nonlinear Poisson-Boltzmann equation for the spheroidal particles. In section IV some 
Monte Carlo simulation results for systems of spheroidal colloids with a point like charge at their centers and explicit 
ions are presented. 
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II. LINEAR POISSON-BOLTZMANN EQUATION 



We consider a single colloidal particle with spheroidal shape, immersed in the electrolyte. We choose the center of 
the coordinate system in the center of the center of the colloidal particle. In spheroidal coordinates, the surface of 
the colloidal particle is given by ^ = ^o- 

In the Debye-Hiickel regime, Poisson-Boltzmann equation for the potential becomes 

V2*(C,r?,0) = /c"*(?,r?,<^) (6) 
= A{£^)B{ri)C{(l)), the equations in each of the 



which is separable in spheroidal coordinates. Writing "i!{(^,r], 
coordinates for the prolate case are 
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with a = na. The solutions for (9) are 



(7) 
(8) 
(9) 

(10) 



The functions of the ^ and 77 coordinates (equations (7) and (8)) satisfy the same differential equation, known as 
the spheroidal equation, which is a generalization of Lcgcndre equation. The only difference between the radial and 
angular functions is the range in which they are defined (— 1 < < 1, and 1 < ^ < 00 in the prolate case, or < ^ < 00 
in the oblate case). The constant of integration A is a function of integers / and to,, and in the spherical (Legendre) 
case it is A = Z(^ + l), for Z = 0, 1, 2, . . .. The solutions (radial and angular) arc expressed as expansions in the Legendre 
functions and are known as the spheroidal wave functions. We use the notation and normalization presented in [15, 16] 
for the spheroidal wave functions. In this work, the numerical calculation of the spheroidal wave functions is done 
with Mathematica, which uses the same conventions as in [15, 16]. The angular functions are denoted ps™. They are 
the analogous to the associate Legendre functions P™ in the spherical case. The normalization used here is [15, 16] 



^ , m/- m2j 2 {l + my. 
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There are two groups of solutions for the radial equation 5™^^^ and S^^^K For 00 these functions behave as the 

spherical Bessel functions 



tat, 



-,m(2) / .- 



(ia,^) ~ ni(m^) ~ — - sin(io^ — i7r/2). 



(12) 
(13) 



To satisfy the boundary condition ^{^, 77, ^) — >■ when ^ — ^ 00, the appropriate radial function is the combination 



(14) 



which is analogous to the spherical Hankel function for the spherical case. The solution to equation (6) is now 
expressed as a superposition of spheroidal wave functions 



l,m 



for the prolate case. In the oblate situation, 

00 

n^, vA) = Yl «re''"Vr(a, »?)i?r («. ^o- 
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From equations (12) and (13), it can be seen that the behavior of the radial functions of equation (14) as ^ ^ oo 
is the same as the one of the screened potential (eq. (1)). Indeed, when ^ — > oo, we have ~ r, and r] ~ cos6, and 

mi~a,0-^- (17) 
Therefore, far from the colloid surface, » 1, the potential behaves as 

„~Kr 

^{r,e,<k)-ziB f{h,e,(k)-^, (18) 

where the anisotropy function /(a, 6, (f)) is given by 

f{a,6) = -^Y.aTe'^'t'psT{i~a,cos6). (19) 
Z Ibk j—' 

In this work, wc will only consider spheroids with a surface charge density which is rotationally invariant around the 
z axis. Therefore the potential will be independent of the azimuthal coordinate (j) and only the terms with m = 
survive in equation (19). 

To characterize the anisotropy along the minor and major axis directions, we define a "maximum anisotropy 

function" 

|/(a,7r/2)-/(5,0)| 
mm[/(o,7r/2),/(a,0)] 

In the next section, wc will (;ompute the anisotropy function for various aspect ratios maintaining the size of the 
major semi-axis (i?) constant. 



A. Constant surface charge density 

A case which allows for an analytical solution to Poisson Boltzmann equation in the linear Debye-Hiickel regime is 
when the colloidal particle is penetrable by the ions of the electrolyte [11] and it has a constant charge density a in 
its surface ^ = ^o- This model can be applied to cells or vacuoles with ion channels, or as a coarse grained description 
of a globule formed by a hydrophobic polyelectrolyte [17] in a poor solvent. The total charge of the colloid is 



(21) 



In this case, the electrostatic potential satisfies the linear Poisson-Boltzmann equation (6), inside and outside the 
spheroid, with the same Debye length k~^. Therefore, the potential inside of the particle has the following form, in 
the prolate case, 

oo 

*™(^,'?) = Y,b^iPs^i{m,n)S°^'\m,(,). (22) 

Here, only the radial function S*™*-^' is used for the radial part, since 5™^^^ diverges as ^ goes to 1 (prolate case) or 
(oblate case). On the other hand, the potential outside of the particle has the form shown in equation (15), with 
m = 0. 

The coefhcients a° and 6° are obtained by imposing the boundary conditions that the electrostatic potential at the 
surface of the spheroid ^ = is continuous, and that the discontinuity of the normal (n) component of the electric 
field at the spheroid surface is proportional to its surface charge density 



From the continuity of the electrostatic potential, we obtain 



A-Klea/e = V*(Co", ??) • n - V*(^^, 77) • n 

(23) 
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Using equation (23) and the orthogonality of the angular spheroidal wave functions 

psi''{ia,r])ps",\ia,f])dr] = ^_^ (| + "^)- j^^, (25) 
.1 it + 1 (t — my. 

we obtain the coefficients 

et » j_i 

The angular spheroidal wave functions have the same parity as therefore for I odd, the coefficient = 0. 
Replacing (26) and (21) into (19), we obtain the anisotropy function for the prolate case 

Ei(2« + l)r'5°(^\m,^o)psO(m,cos^)/^,Vir^P«?(i«>^)rf^ 
/(a, 0)pro = 772 T^c2, -1/1/ • 

In the oblate case, it is 

f{a, e)ou = ; T-7fT^ • (2^) 

The sum over / runs through all even integers. 

Figure 2 shows the anisotropy function for the constant surface charge case as the aspect ratio (R'/R) changes from 
the sphere to the rod for the prolate case and from the sphere to the disc for the oblate case. 

For a prolate spheroid, the anisotropy function, and thus the electrostatic potential at large distances, is larger in 
the direction ^ = of the large axis. For an oblate spheroid, the anisotropy function is also larger in the direction of 
the large axis (now labeled 6 = 7r/2). 




FIG. 2: Anisotropy function for the oblate(top) and prolate(bottom) cases for various aspect ratios {R'/R) and a constant 
surface charge density. Here the major semi-axis R is five times the Debye length. 

Now, we study the difference in the potential along the major and minor axes directions (^ = and 9 = 7r/2), 
with the maximum anisotropy function defined in Eq. (20). The maximum anisotropy function for various values of 
the Debye length and different aspect ratios is shown on figure 3. It can be seen that the function is small when the 
Debye length is greater than the dimensions of the particle, as the big cloud of ions screens the charge on the colloid 
surface in such a way that far from the colloid the potential is almost isotropic. On the other hand, for small values 
of the Debye length the potential retains its anisotropy even at large distances. It can also be seen that the shape of 
/m vs. aspect ratio changes its curvature as the Debye length is increased. 
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FIG. 3: Mciximum anisotropy function for the oblate and prolate cases for various aspect ratios and Debye length values with 
a constant surface charge density. 

Another boundary condition of practical interest is when the colloidal particle has constant surface charge density 
but it is impenetrable by the electrolyte. Unfortunately, this situation turns out to be more difficult to treat analyti- 
cally. The reason is the following. Inside the spheroidal particle, the electrostatic potential satisfies Laplace equation. 
In spheroidal coordinates, the corresponding solution to Laplace equation is [14] 

*in(e,r/)=5]diPi(emW (29) 

I 

where Pi are the Legendre polynomials. The potential inside (29) and outside (15) are expanded in different basis of 
functions for the angular coordinate rj. To satisfy the boundary conditions at the surface ^ = of the spheroid and 
find the coefficients di and , one has first to change the basis of angular functions to use the same basis inside and 
outside, i.e. express the Legendre polynomials in terms of the spheroidal wave functions ps'i or vice-versa. We will 
not attempt to follow this route. For most practical situations, the colloidal particle will be highly charged. In that 
situation, the nonlinear Poisson-Boltzmann equation should be used. As we will explain in the following section and 
in section III, for highly charged colloids, the linear Poisson-Boltzmann equation with Dirichlet boundary conditions 
(fixed constant potential at the surface of the particle) is more relevant. 

B. Constant surface potential 

Now, we turn our attention to the case when the colloidal particle is maintained at a fixed constant potential and 
it is not penetrable by the electrolyte. This situation is relevant for the case of highly charged colloidal particles. 
Indeed, as explained in [8, 18, 19], for highly charged colloidal particles, the constant potential boundary condition 
can be seen as an effective boundary condition to join the solution of the linear Poisson-Boltzmann equation valid 
far from the colloid, with the solution of the nonlinear Poisson-Boltzmann equation valid near the colloid surface. 
As a first approximation, the value ^'o for the surface potential is obtained using the known solution of nonlinear 
Poisson-Boltzmann equation for a planar highly charged surface. For a charge-symmetric electrolyte it is Vfo = 4. 

Outside the colloidal particle, the potential is given, for a prolate spheroid, by equation (15), or, for an oblate 
spheroid, by equation (16). Imposing ^'(Co,'?) — ^'o, and using the orthogonality of the angular spheroidal wave 
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functions, we obtain the coefficients 



21 + 1 



*o^7Do77m / (30) 



in the prolate case. Due to the parity properties of ps", if I is odd, then aj* = 0. 
The anisotropic surface charge density of the colloidal particle is 



^^a°<(m,^oK(m,77), (31) 



47rZBa ^eFV 

where the prime denotes the derivative with respect to ^. Integrating over the surface of the colloidal particle, we 
obtain its total charge 

ZIb = -^^^^=^^a°i?o'(m,^o) f ^ps°i{i~a,7i)drj, (32) 



Let us define 



Ci{z) = J ^ps'^{z,v)dv, (33) 



which is nonzero only when / is even. Then, the charge; of the colloid is 

aM^i - 1) ^ 1 ^r< ^,;;^2^RL^!Mo) 



ZIb = — 



Finally, the anisotropy function, for the prolate situation is 



^ ~'~ -Ci{ia) ps'i{ia, cos 6) 



2 ^R?{ia,^o) 

f{a,9),ro - -5(^2 _i) ^(2/ + i)Q(za)^i?r(-,Co) ' ^''^ 



i?0(w,eo) 



and, for the oblate case, it is 



2 

- "a(gg + l) ^ (2/ + l)q(.)^<(»./^o) • ^''^ 

Figure 4 shows the anisotropy function for different aspect ratios R'/R for the constant surface potential case. 

Again, the anisotropy function is larger in the direction of the major axis of the spheroid, as it was in the case of the 
constant surface charge density. However, the overall magnitude of the anisotropy function is larger in the present 
case of constant surface potential than in the constant surface charge situation of the previous section. This can be 
understood as a "tip" effect. Indeed, with a constant surface potential boundary condition, the tips of the spheroids 
acquire a larger surface charge density, which translate into a higher anisotropy function in that direction (the large 
axis direction). 

Figure 5 shows the maximum anisotropy function for various values of the Debye length and different aspect ratios. 
Comparing it to figure 3 it can be seen that for the present case when the potential at the surface of the colloid is 
constant the values of the maximum anisotropy function for the prolate and oblate cases become different for small 
Debye lengths (compared to the large semi-axis) and for highly anisotropic objects, i.e. small R'/R values, that is, 
objects near the limiting cases of the disc or the rod. 

Figure 6 compares the maximum anisotropy for the oblate and prolate cases for a large value of kR, and for the 
two different situations of constant surface charge density or constant surface potential. It is clear that at constant 
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FIG. 4: Anisotropy function for the oblate(top) and prolate(bottom) cases for various aspect ratios (R' /R) and a constant 
surface potential. Here the major semi-axis is five times the Debye length. 
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FIG. 5: Maximum anisotropy function for the oblate and prolate cases for various aspect ratios and Debye length values with 
a constant potential at the surface. 

surface charge density the function Jm behaves similarly for the prolate and oblate cases while at constant surface 
potential the anisotropy is greater in the prolate case. This is probably related again to a "tip" effect, as mentioned 
earlier. In the Dirichlet boundary condition case, the surface charge at the tips will be more important for prolate 
spheroids than for oblate spheroids, since the local curvature at the tip is higher in the prolate case. This translate 
into a higher anisotropy function in the large axis direction for prolate spheroids, as shown in figure 6. On the other 
hand, the constant surface charge density situation does not discriminate between prolate and oblate spheroids as far 
as the anisotropy function is concerned. 
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FIG. 6: Mciximum anisotropy function for the oblate and prolate cases for kR = 10. 



C. Special cases: the disc and the rod 



In this section, we give some additional details on two special cases, the disc and the rod geometries. A disc can 
be seen as an oblate spheroid with its surface given by = 0. Disc shaped colloidal particles have been studied 
in [8-10, 20, 21] as they are an appropriate model for clay suspensions, like laponite. In [10], the screened electrostatic 
potential was studied within the linear Poisson-Boltzmann theory, with Neumann boundary conditions, i.e. constant 
surface charge density on the disc. At large distances from the disc, compared to the Debye length, the potential is 
of the form (18): a Yukawa screened potential, multiplied by an anisotropy function. The anisotropy function was 
computed exactly in [10], for this situation 

fl~ a\ - 2-^1 (a sing) 

/ l^a, Cjdisc. Neumann — r — i — 1\ ■ l.'J ' J 

a sm 

Here a — ko, with a the radius of the disc, and Ii the modified Bessel function of order 1. Comparing to the general 
expression (28) for oblate spheroids, and taking the limit 0^, we obtain an interesting expansion formula for the 
modified Bessel function in terms of spheroidal wave functions 



£ {2l + l)\f W\ps°i{h,ri')d^' 

-n 9 L-'O 



1=0,2,. 



i-'S'i''\iQ+, a)psna, v) = ^T^r—^ ' ■ (38) 
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For the case of a disc held at constant surface potential, our work gives additional information to the one found in 
the current literature [8]. In [8], this problem was studied, using cylindrical coordinates. Unfortunately, in cylindrical 
coordinates, the linear Poisson-Boltzmann equation turns out to be untractable analytically, when the potential at 
the surface of the disc is fixed. The technical difficulty is that the problem involves mixed boundary conditions 
when cylindrical coordinates (p, 0, z) are used: a Dirichlet boundary condition for p < a and z = 0, and a Neumann 
boundary condition for p > a and z = 0. Using oblate spheroidal coordinates, this problem disappears, only a 
Dirichlet boundary condition is imposed at ^ = 0. Specializing equation (36) for the disc case (^o ~^ O"*^); we obtain 
an analytic expression for the anisotropy function for the disc shaped colloidal particles at fixed potential 

/(a, 0)disc, Dirichlet — 7ol~rTTF<27-Tiri307^ ■n+\ ' (39) 

a s—^ [21 + l)C[{a) o^R\ [a, lO^) 
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An analytical expression for the anisotropy function in this situation was previously unknown. Exactly, it was only 
known that [8] f{a,0)iiisc,Dirichiet = 1- This can be verified numerically in the previous expression, and gives yet 
another interesting relation between the spheroidal wave functions 



ps^(l,a) = -- J ^ps^{d,r])dr] ^^{d,i^) 



(40) 



Since an analytical expression for the anisotropy function in this case was previously unavailable, the authors of [8] 
proposed an approximate analytical ansatz based on a two parameter fit, equation (16) of Rof. [8]. In figure 7, we 
compare the ansatz from [8], with the exact analytical expression (39). As it can be seen in the figure, the two 
parameter ansatz proposed in [8], is a fairly good approximation. It underestimates the anisotropy, but, for a = 5, 
the relative maximum error is for 6 = 7r/2, and it is only 5%. 




FIG. 7: Anisotropy function for a disc at constant surface potential with radius 5 times the Debye length. 



The screening of a rod- like stiff polyelectrolyte, such DNA fragments. Tobacco Mosaic Virus (TVM), etc., can be 
described by our general approach, since a rod can be seen as a prolate spheroid with — 1. In [13, 22], the screened 
electrostatic potential and the anisotropy function of rod-like polyelectrolytes, modeled as cylinders of finite length 
and nonzero radius, were computed numerically. From our approach, we obtain analytical results for the limiting case 
of infinitely thin rods (radius zero). For a thin rod with fixed charge density, the anisotropy function can be obtained 
from (27), taking the limit ^ 1, 



f{a,0)r 



)d,N. 



cumann 



= - V(2; + i)z- 



■'5r«(l,za) 



Vl-7?'V°(a,r;')d7?' 



ps°(a, cos( 



(41) 



Here 2a is the length of the rod. 

The case of a rod at fixed constant potential is somehow pathological, as the spheroidal wave function i?['(za, ^o) 
diverges as ~^ 1- This difficulty already happens in the absence of the electrolyte. A charged rod in the vacuum 
having a fixed constant potential would have an infinite linear charge density. Indeed, in the vacuum, the potential 
of a prolate spheroid with fixed potential on its surface (^ = ^o) is [14] 



^prolate, vacuum (0 — ^0 



In^ 

&-1 



This expression has a logarithmic divergence when ^ 1- The corresponding surface charge density is 



0-(??)v 



(42) 



(43) 



In the limit of a rod, 1; one can define 

surface element is dS = a}\J (^g ~ ^^)(€o ~ 1) 



the corresponding linear charge density A(r?), remembering that the 



= lim ^aJ{eo-v'){eo-^)<^{v) = Hm 



'I'n 



ln|2+i 

5o-l 



(44) 
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The linear charge density has the same logarithmic divergence when ^ — > as the potential. 

This same difficulty appears in the presence of an electrolyte for the screened potential, the potential and the charge 
have a logarithmic divergence when 1- Nevertheless, the anisotropy function will have a finite limit when .^o 1, 
since it is defined as the ratio between the potential at large distance and the total charge of the rod: the logarithmic 
divergences will cancel. When — ^ 1) 

i??(^a,6)~5/(a)ln(eo-l). (45) 
We shall determine the explicit expression for gi{a) later (equation (48)). The derivative with respect to ^ diverges as 

<(ia,^o)-|^ (46) 

Therefore, in the limit — ^ 1) the anisotropy function (35) has a finite limit 

ft- j:i{2l + l)Ci{z~a)ps^ita,cos9)/gi{a) 

f{a, ^).od,Diochiet YM+T)Cfm ^ ' 

with Ciiia) given by (33). The prefactor giia) in the asymptotic formula (45) can be determined as follows. In the 
direction perpendicular to the rod, 6 = 7r/2, the anisotropy function is equal to 1. Therefore 



gi{a) = - .'J; . 48 
a Ci {la) 



Finally, 



/ (a, 0).od,DiHchiet = J2^(2l + l)Q{dr ' ■ ^^^^ 

As a subproduct of this analysis, we obtained the asymptotic behavior of the spheroidal wave function, when ^ — >■ 1, 

~ -f^ln(e-l). (50) 



III. NONLINEAR POISSON-BOLTZMANN EQUATION 

A. Renormalized anisotropy function 

When the colloidal particle is highly charged, Poisson-Boltzmann equation cannot be linearized. The full nonlinear 
equation reads 

V^* = [e^-* - e-^+*] . (51) 

z+ + Z- 

Although near the surface of the highly charged colloidal particle, the full nonlinear equation should be used, at 
large distances from it, the potential becomes small due to the screening. Therefore, at large distances, compared 
to the Debye length, from the colloidal particle, the potential satisfies again the linear equation (6). For a spherical 
symmetric colloidal particle, this means that the potential at large distances from the particle will be of the DLVO 
form (1) 

*W~^.en/B(^^^, (52) 

but now, Zren IS uot the bare charge Z of the colloidal particle. This new prefactor to the DLVO potential is known 
as the renormalized charge [18, 19, 23]. 

For spheroidal colloidal particles, we expect that, at large distances, the potential behaves as (18) 

p-nr 

^{r,e,^)^F{h,e,^) . (53) 
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The determination the anisotropy function F requires to apply the boundary condition at the surface of the coUoid. 
But in this region the hnear solution is not valid. Therefore the anisotropy function will be in principle different 
from the one obtained in the previous section. For a spherical symmetric colloidal particle, the anisotropy function 
is constant, and in the linear case, for impenetrable spheres of radius R it is /o = e'^^/{l + kR). Writing down the 
anisotropy function F (here constant) for the nonlinear situation as 

F = Z,,jB.fo , (54) 

defines a renormalized charge Z^-en for the colloidal particle, which encodes all the nonlinear phenomena due to its 
large charge. In the anisotropic case, one might be tempted to write 

F{a, e, (t>) = Z^^JeMa, 9, </>) . (55) 

However, there is no guaranty that the nonlinear effects would only affect the anisotropy function by an overall 
multiplicative factor (the renormalized charge), they would probably also affect the dependence on the angular vari- 
ables {0,(j)), i.e. fni{a,9,(j)) will be different from its linear counterpart. Therefore the full anisotropy function F is 
"renormalized" by the nonlinear effects, and the separation proposed in (55) is not unique. 

For disc shaped colloidal particles, in [8], the following approach was used. Since the linear anisotropy function 
satisfies /(a,0) = 1 in the direction perpendicular to the disc, the same relation is proposed to define nonlinear 
anisotropy function /ni(a,0) = 1. From this, a unique definition of the renormalized charge follows, using Eq. (55) 
when 9 = 0. From a numerical resolution of the nonlinear Poisson-Boltzmann equation one can obtain the anisotropy 
function for the disc in the nonlinear regime. 

For spheroids other that the disc or the rod, the property /(a, 6) = 1 does not hold for any direction 6. Therefore 
there is not an unambiguous way to define a renormalized charge Z^en and a nonlinear anisotropic function /ni . Instead 
of trying to separate the nonlinear effects into a renormalized charge and an anisotropy function, we will study the 
combined effect in the complete anisotropy function F defined by equation (53). 

To find the explicit form of the anisotropy function F, in principle, one needs to solve the nonlinear equation 
and match the linear solution (53) valid only at large distances, with the nonlinear one. This can always be done 
numerically. However an analytical approach was proposed in [8, 18, 19], valid for highly charged colloidal particles 
with dimensions larger than the Debye length {kR ^ 1). Under these conditions, the region where the nonlinear 
equation should be used is very close to the colloidal particle, and at such close distances, the colloidal particle can 
be seen as a planar object (as a first approximation) . Therefore one can approximate the nonlinear potential in that 
region with the known planar solution of the nonlinear Poisson-Boltzmann equation in one dimension. Then, the 
linear potential valid at larger distances is matched to the nonlinear one through an effective boundary condition of 
the Dirichlet type at the surface of the colloidal particle, i.e. ^"(^0, 0) = ^'o- The effective value \I'o of the potential 
is obtained from the solution of the nonlinear one-dimensional Poisson-Boltzmann equation. Provided that the charge 
of colloidal particle is large, ^'o does not depend on the charge of the colloidal particle, only on the constitution of 
the electrolyte, i.e. on the valencies z+ and Z-. For a charge-symmetric electrolyte, 2+ = z_ = 1, it is \I/o — 4. For 
further details on this picture of highly charged colloidal particles seen as objects at constant potential, we refer the 
reader to references [8, 18, 19]. 

In section II B, we solved the linear Poisson-Boltzmann equation with the constant surface potential boundary 
condition. From the analysis of that section, we obtain directly an analytic expression for the renormalized anisotropy 
function. For prolate spheroids, 

21 + 1 

F{d,9)pro = ^oV „ „o/~ ^ . Ci{ia)ps^ {ia, COS 9), (56) 
^ 2KR\'(ia,^o) 

and for oblate spheroids, 

F{a,0Ui = ojlt Ciia)ps1 {a, cos 9). (57) 

In the following section, we test the constant surface potential boundary condition, by solving numerically the 
nonlinear Poisson-Boltzmann equation for spheroids. 

B. Numerical resolution of Poisson-Boltzmann equation 

We consider a single spheroid with surface charge density a immersed in an electrolyte, and the spheroid is penetrable 

by the electrolyte. In the mean field approximation, the screened electrostatic potential satisfies the nonlinear Poisson- 
Boltzmann equation (51) outside and inside the spheroid. This differential equation is complemented by the boundary 
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conditions V^'(^, 77, </>) when ^ — > 00, ^' continuous at the surface of the spheroid £, = (.0, and the discontinuity of 
the electric field at the surface of the spheroid is proportional to cr, Eq. (23). To solve equation (51) numerically, we 
implemented an algorithm similar to the one used in Refs. [8, 22] for rods and discs. 

Let Pcoi(C) = "■'^(C ~ Co)/^c be the charge density of the colloidal particle, with the scale factor corresponding 
to the ^ coordinate [14], = a^J (^^ — rf')/{£,'^ — 1) for prolate spheroids, /ij = a^J {(^"^ +'7^)/(^^ + 1) for oblate 
spheroids. Poisson-Boltzmann equation (51), with the boundary conditions imposed on ^, can be rewritten as 



k2 



V^* - k''^ = [e-* - e-^+*] - - MlB/e)p.oi ■ (58) 

Z-^- + Z— 



Introducing the kernel G of V — k in free space 



1 „-re|r-r'| 

G(r,r') = - — ^ ^, (59) 

47r r — r 



equation (58) can be cast as an integral equation 



^-(r)= / G(r,r' 

iM3 



Z+ 



gZ_*(r') _ g-^+*(r') 



- K^^{r') - 47r(/s/e)pcoi(r') 



dr' (60) 



The numerical algorithm to solve this equation is iterative. It starts with a trial value for the potential, computes 
the integral in the right-hand side of equation (60) to obtain a new value for the potential, ^'out; then uses a linear 
mixing of the two, (1 — a) ^'in + agouti to start the next iteration, with the parameter a in the range 10~^ - 10~^. This 
procedure is iterated until convergence is achieved, i.e., / |^'in — ^'outl/ / |^in| < e, with e a relative error tolerance of 
the order 10^'''. 

For the present problem, we used spheroidal coordinates. The kernel G can be expanded in spheroidal wave functions 

as 

G(e,r,,</.;e',r?',0') = -^J2^2l + l)e"^^^'f'''f''^psY'il~a,v)ps^{^&,v')^''S^^'\ (61) 

l,m 

for a prolate spheroid, and 

G(C, rj, 0; e, V', 0') = "^^ E(2^ + l)e^'"'^"^'Vsr(a, v)psr{a, sf'\a, iU)RTifl. *^>) (62) 

for an oblate spheroid, with ^< = min(^,^') and ^> = max(^, ^'). Note that, since the charge density of the colloidal 
particle does not depend on the azimuthal angle, in the integral equation (60), only the terms for m = remain after 
the integral over the azimuthal angle </>' is performed. 

Numerically, one cannot compute the integral in equation (60) in the whole space M^, but rather it has to be 
restricted to a finite region ^ < ^max- The integral equation (60) with the integral restricted to a finite region ^ < ^max, 
is equivalent to the nonlinear Poisson-Boltzmann differential equation (58), for ^ < ^maxj with the boundary condition 
that outside the region of integration (^ > ^max) the potential satisfies the linear Poisson-Boltzmann equation (6). 
This is perfectly acceptable provided that a^max ^ 1, since due to the screening we know that the potential will be 
small, \1/ ^ 1, at large distances from the colloidal particle, and thus it will satisfy the linear Poisson-Boltzmann 
equation. 



C. Numerical results 



Figure 8 shows the numeric solution of the nonlinear Poisson-Boltzmann equation for a prolate spheroid with 
^0 — 1-2 and a — 8.0, kR = 9.6, kR' = 5.3, and aspect ratio R'/R = 0.55. Its surface charge density is given by 
lBcr/{K,e) = 10. The corresponding bare charge is then given by Zls/a o± 690 :» 1, which is large, thus in the nonlinear 
and saturated regime. The spheroid is immersed in a charge-symmetric electrolyte with z+ = Z- = 1. In the same 
figure, we show the solution of the linear Poisson-Boltzmann equation, with the effective boundary condition of fixed 
potential 4*0 = 4 at the surface of the spheroid 

*iin(C,»7) = ^o^ ^^o,t\ . Ci{ia)psf{i~a,Tj)Rnia,0 ■ (63) 
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FIG. 8: Electrostatic potential for a prolate spheroid with = 1-2, a = 8.0, and surface charge density given by lB<y/{Ke) = 10. 
The symbols are from a numerical resolution of the nonlinear Poisson-Boltzmann equation, and the lines from the solution of 
the linear Poisson-Boltzmann equation with the fixed surface potential boundary condition. 



At larges distances from the spheroid, Ka^ > 2, the linear approximation is excellent. This strongly supports the 
picture of constant potential objects for highly charged colloidal particles presented in the previous section and in 
Refs. [8, 18, 19]. 

In figure 9, wc present the results for an oblate spheroid, now with = 0.50, a = 3.0, kR' = 1.5, kR = 3.4, 
R'/R = 0.45 and charge density lBcr/{K,e) = 10. The agreement with the linear solution for fixed surface potential '^q 
is again very good. Since the dimensions of the spheroid are not very large compared to the Debye length, one could 
eventually search for size dependent (in kR and kR') corrections to the planar value ^'o = 4 to improve the linear 
approximation. As an additional test, in figure 10 and 11, we present the results for a charge-asymmetric electrolytes 
with z+ = 2, Z- = 1, (2:1), and the reverse situation Zj^ = 1, Z- = 2, (1:2). In these cases, the saturation values 
of the effective potential are = 6 for the (2:1) situation, and = 6(2 — for the (1:2) situation [24]. Again, 
the agreement between the nonlinear solution and the linear one far from the colloidal particle is excellent, in all 
directions (in the figures we report only the results for the directions t] = 0, i.e. 9 = ■jt/2 and t] = 1, i.e. 6 = 0). 




FIG. 9: Electrostatic potential for an oblate spheroid with ^0 = 0.50, 2 = 3.0, and surface charge density given by lBO'/{ne) = 10. 
The electrolyte is charge-symmetric with valencies z+ = 1 and Z- = 1. The symbols are from a numerical resolution of the 
nonlinear Poisson-Boltzmann equation, and the lines from the solution of the linear Poisson-Boltzmann equation with the 
fixed surface potential boundary condition. 
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FIG. 10: Electrostatic potential for an oblate spheroid with = 0.50, a = 3.0, and surface charge density given by lBcr/{Ke) = 
10. The electrolyte has valencies z+ = 2 and Z- = 1. The symbols arc from a numerical resolution of the nonlinear Poisson- 
Boltzmann equation, and the lines from the solution of the linear Poisson-Boltzmann equation with the fixed surface potential 
boundary condition. 




FIG. 11: Electrostatic potential for an oblate spheroid with = 0.50, a — 3.0, and surface charge density given by /_Ba/(Ke) = 
10. The electrolyte has valencies z+ = 1 and Z- = 2. The symbols are from a numerical resolution of the nonlinear Poisson- 
Boltzmann equation, and the lines from the solution of the linear Poisson-Boltzmann equation with the fixed surface potential 
boundary condition. 



IV. MONTE CARLO SIMULATIONS 



In an attempt to study the charged spheroidal colloid system beyond mean field theories, we made some preliminary 
simulations of a system of hard spheroidal colloids with a point charge Ze at their centers, and explicit colons and 
counter-ions of charge +e and — e confined in a simulation box of side length L, and using periodic boundary conditions. 
The reason to choose the model with a charge at the center of the particles rather than with constant potential or 
charge density at their surface was a matter of computational efficiency, as in the latter cases the analytic expression 
for the; interaction between the particles is very complicated (except for the spherical case) and thus the potential 
needs to be computed numerically. On top of that, as we use periodic boundary conditions, it would be necessary to 
calculate the Ewald sums form of the potential to take correctly into account the images of the system, which makes 
it even more difficult, and numerically time consuming. 

Although the simulations were carried out at relatively low colloid densities, with a colloid volume fraction of 
0.0025, the systems are not at the infinite dilution limit, and five colloids were simulated in the simulation box, rather 
than one, to avoid the bias that a crystalline structure could introduce to the system. Unfortunately, the analytical 
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solution of Poisson-Boltzmann equation in the case of point charges at the center of the spheroidal colloids turns out 
to be difficult even in its linearized (Debye-Hiickel) form, as the symmetry of the potential inside and outside of the 
spheroidal particle is not the same. For these reasons, we will not attempt a direct quantitative comparison of the 
simulation results with the solution of Poisson-Boltzmann equation, but rather to make a qualitative analysis. 



A. Model 



The particles interact at short range with the spheroidal hard core potential, and at long range with the Coulomb 
electrostatic potential. The dimensionless long range Coulomb interaction potential {(3ex potential) between particles 
i and j is given by 



Uc{rij) 



ZiZjlB 



(64) 



where Zi and zj are the valences of particles i and j and rij = \rj — ri\ is the distance between the centers of the 

particles. As we use periodic boundary conditions, in order to take correctly into account the interactions between 
particles in the main simulation box and particles in the periodic images, Ewald sums were used, so the Coulomb 
potential energy takes the form [25] 



Ec 



ericjanj) 

ZiZjtB Z 



1 

Z3 



(65) 



where n = {nx,ny,nz) G Z"^ is the image vector in real space, fj 



Ln\, the prime denotes that i ^ j if n = 0, 



k = 27rp/L, with p e Z^, is the grid vector in the Fourier space and a is a parameter controlling the convergence of 
the real and Fourier space sums. We use the minimum image convention so that n = 0. The Fourier sum is truncated 

at IpI = 9 with a = 0.76. 



The short ranged (hard-core) potential is 

l^Hc{^ij,^i, ^i) 



0; F(r,j,f7„%) > 1 
oo; F(rij,0i,0j) < 1 



(66) 



where rjj is the separation vector between the centers of the ellipsoids and F{rij,Qi,Qj) is the contact function 
between two ellipsoids with orientations Oj and Qj as defined in [26]. This function F has the property that for two 
ellipsoids A and B 



F{rAB,^A,^B) 



< 1 if A and B overlap 

= 1 if A and B are externally tangent 

> 1 if A and B do not overlap 



(67) 



Each spheroidal particle has a point charge at its center. Alternatively, one can think that each spheroid is void 
inside and has a certain charge density in its surface that creates an electrostatic potential outside the spheroid, at 
a distance r from its center, proportional to 1/r. Now, we compute this equivalent surface charge, which could be 
used in principle as a Neumann boundary condition to solve analytically the Poisson-Boltzmann equation. The Green 
function of Laplace equation, in spheroidal coordinates [14] , is 



1 



^ OO n 



n=0 



m=0 



(n — m)l 
(n -|- m)! 



cos[m{(j) ■ 



xp™(,7')^r(^) 



(68) 



in the prolate case, and 



m+l 



-E(2-+i)E 

m=0 

xP™(V)P„"(r?) 



n=0 



COS ml 



(n — m)! 
(n -|- m)! 

/ p™(ic')Q;r(*e); e > c 



(69) 
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in the oblate case, with P„ and Qn the Legendre functions of degree n of the first and second kind respectively. The 
value of is 1 if TO = 0, and 2 otherwise. For the case of spheroids of revolution, and taking r' at the origin, we 
have that only the m = terms are non zero. Then, the potential inside and outside of a spheroidal particle is 

Uout{^,v) = ^ = ^ V(2n+ l)P„(0)P„(r?)Q„(O, (70) 
r a ■^—i 

oc 

Kin{^,V) = yAnPn{v)Pn{0- (71) 



n=0 



in the prolate case, and 



Uouti^,v) = ^ = V(2n+ l)P„(0)P„(»7)Q„(iO. (72) 

OO 

{v)Pn{iO- (73) 

n=0 

in the oblate case. The constants are determined by using the boundary condition 

Kin{^o)=Uout{Co)- (74) 

Then, using 



1 dU, 



out 



1 dUi, 



-4./.^ (75) 



€=io 



along with equations (70) to (73), it is possible to find a surface charge density equivalent to having a point charge Ze 
at the center of the spheroids (i.e. the charge density that produces a potential outside the spheroidal particle that is 
proportional to 1/r) 



a{r]) 

for the prolate case, and 



Ze 1 



V^a«^iV ^^n (2n-l)!! P2n{v) 



^^^^ VW^WWl^) fo (2n)!! P2„(?o) 



^'^~^--'VWTWmTT)^r^'^^ (2n)!! P.„(^6)• ^''^ 

for the oblate case. Here the difficulty of solving analytically the linear PB equation (eq. (6)) becomes evident, as 
the surface charge density is given in terms of the Legendre functions, while solution of the Helmholtz equation in 
spheroidal coordinates are the spheroidal wave functions. It is worth mentioning here that this charge distribution, 
which produces a radial field outside of the particle and the charge distribution associated with a constant surface 
potential (eq. (31)) have different behaviors (Figure 12): the radial field charge distribution is higher in the direction 
of the small axis, while constant potential surface charge is higher in the direction of the large axis. 



B. Simulations 



Monte Carlo simulations were carried out for prolate and oblate colloids with R' / R = 1/3, and three different 
values of = —50, —200 and —400 for each case. The Bjerrum length used is = 0.0625 (in this section the lengths 
are normalized as r* = r/{8RR'^y^^ in the prolate case and r* = r/{8R'R'^y^^ in the oblate case). Each system 
consisted of 5 spheroidal colloidal particles, with a packing fraction rj = 0.0025, 500 colons and 750, 1500 and 2500 
counterions respectively. The reduced volume is V* = 1028.19. 

From the simulations the colloid-counterion correlation functions h = g — 1 

h{nj, Qi, Qj) = E E E /^'""'(r-«)^So"'(fii. ^jM (78) 

m n I 
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FIG. 12: Surface charge density as a function of rj for the oblate {R'/R = 1/3, upper) and prolate [R' / R = 1/3, lower) cases. 
The constant potential case was calculated with ^'o = 4 and the radial field case with Ze — 30. 



were obtained. The interparticle distance is r^j, fli are the Euler angles that give the orientation of particle i, R is 
the orientation of the interparticle vector, and the are the rotational invariants defined in [27] multiplied by a 

prefactor 

r"' = ;!/(^™ "^[y (79) 

where the term on the denominator is a Wigner's 3j symbol. 

In our case, since the microions have spherical symmetry and the spheroids have rotational symmetry around their 
revolution axis as well as the perpendicular plane to this axis that crosses at their center, the expression (78) reduces 
to a simple expansion in Legendre polynomials 

oo 

hir,e) = Y,h^'ir)P2i{cos0), (80) 

where 9 is the angle of the interparticle vector with respect to the orientation of the revolution axis of the colloid and 

/i2'(r) = l^ili /" P2i{cos9)h{r,e)d{cos9). (81) 

The radial functions h?^{r) were obtained from the simulations and then used to reconstruct the correlation function 
using equation (80), truncating the summation at Imax = 100. 

The potential of mean force, Umf, associated with the colloid plus its screening cloud of ions was obtained from the 
correlation functions using 

/i(r,0) =e-"'"^('''^^-l. (82) 

The simulations with 1250-3000 ions were equilibrated for 100000-150000 Monte Carlo steps and averaged for 60000- 
130000 Monte Carlo steps. 

The potential of mean force is plotted against r* in figures 13, 14 and 15. In the figures it can be seen that for 
the oblate case the potential increases as we go from the poles of the colloidal particles {9 = 0) to the equatorial line 
{9 = 7r/2) for a fixed value of r* . For the weak coupling case {Z = 50) the potential depends little on the orientation, 
as could be expected for the isotropic (spherical) case, but the effect is more evident as Z increases. This is the same 
behavior that was observed for the solution of the linearized Poisson-Boltzmann in the constant charge and constant 
potential cases (figures 2 and 4), even though the distribution of charge is very different in all cases (figure 12). 
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FIG. 13: Effective potential vs. r for the oblate (up) and prolate (down) cases with R' /R = 1/3 for different angles with 
Z = 50. 
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FIG. 14: Effective potential vs. r for the oblate (up) and prolate (down) cases with R! /R =1/3 for different angles with 
Z = 200. 

In the prolate case the potential decreases as we go from = to 6' = 7r/2 for fixed r* , and it is again the same 
behavior observed for the solution of the linearized Poisson-Boltzmann in the constant charge and constant potential 
cases. 

From these results it can be seen that the short range steric interaction (the shape of the excluded volume) has a 
strong effect on the effective screened potential of the colloid, producing an anisotropy which is larger in the orientation 
where the curvature of the colloid is higher, even for different surface charge distributions. 

V. CONCLUSION 

We studied the effective electrostatic potential created by a spheroidal colloidal particle and its screening cloud, 
within the mean field approximation, using Poisson Botzmann equation in its linear and nonlinear forms. Also, we 
did a preliminary exploration beyond the mean field with Monte Carlo simulations. 

In all the studied cases, we confirmed the persistence of anisotropy of the effective potential at large distances, as 
it was already noticed for other anisotropic particles such as cylinders and discs [8, 13]. For highly charged colloidal 
particles, but within the validity domain of mean field, we wore able to test the picture of constant potential objects 
for these anisotropic particles [18, 19]: at large distances from the charged particle, the nonlinear solution of Poisson- 
Boltzmann equation can be approximated by the linear one with an effective boundary condition of constant potential 
at the surface of the particle. In this anisotropic situation, our work has provided a strong test for this picture, since 
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FIG. 15: Effective potential vs. r for the oblate (up) and prolate (down) cases with R' /R =1/3 for different angles with 
Z = 400. 



a constant potential boundary condition is very different from a constant surface charge density, contrary to the 
situation in the spherical case. 

With this constant potential boundary condition, we found that, at large distances from the spheroidal particle, 
the potential is larger in the direction where the curvature of the particle is higher, that is the large axis direction. 
From the simulations, we observed that the potential of mean force around the spheroidal colloidal particles with a 
point charge a their centers, has a similar behavior to that of the constant surface charge density and constant surface 
potential cases computed analytically, even though the equivalent surface charge density is very different from the 
two previous cases. 
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